%% This script was originally sent by Alexei Kuvshinov
% example of multi-taper analysis

ncid = netcdf.open('/data/backup/mnair/obs_mag_data/update_2011/secular_removed/GNA_1995_2010.nc','NOWRITE');
a = (datenum(2011,1,1)- datenum(1995,1,1))*24*60;

time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2010,12,31,23,59,30))';

X_ID = netcdf.inqVarID(ncid,'Magnetic_Field_X');
Y_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Y');
Z_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Z');

z_data = netcdf.getVar(ncid,Z_ID,0,a);

z_data = double(z_data)/10;
z_data(z_data == 99999.9 ) = NaN;
L = isnan(z_data);
z_data(L) = 0;

[P, F] = pmtm(z_data, 4, length(z_data), 1/60.);






%periods = 3:0.1:100;
periods = [12.421,24];
spectra = [];
for i = 1:length(periods),

x = [cos(2*pi*time_array_min/( periods(i)/24)) sin(2*pi*time_array_min/( periods(i)/24))];	 
 
spectra_model_rem(i,:) = robustfit(x, z_data);

end;